share <- function(delta, mui){
  num   <- exp(delta + mui) 
  den   <- t(sapply(unique(t), function(m) colSums(num[which(t == m), ])))[mkti, ] 
  si    <- num / (1 + den)  
  s     <- rowMeans(si)     
  return(list(si = si, s = s))
}

share_2 <- function(p, deltar, v, zeta){
  alpha <- zeta[2]
  mui   <- - zeta[1]*v*p
  delta <- deltar - alpha*p
  s <- share(delta = delta, mui = mui)
  si  <- s$si 
  s   <- s$s
  return(s)
}

fixedn <- function(delta, mui) {
  s <- share(delta = delta, mui = mui)  
  s <- s$s
  delta <- as.numeric(delta + log(S) - log(s))
  return(delta)
}

fixed <- function(zeta) {
  mui   <- zeta[1]*v*p  
  delta <- try(squarem(delta0, fixedn, mui = mui, control = list(tol = 1e-10)), silent = T)
  expl  <- class(delta)

  if(expl == 'try-error') {
    si <- delta <- NULL 
  } else {
    delta <- delta$par 
    si <- share(delta = delta, mui = mui)
    si <- si$si
  }
  return(list(delta = delta, si = si, expl = expl))
}

struct <- function(zeta) {
  Fout  <- fixed(zeta)
  expl  <- Fout$expl
  
  if(expl == "try-error") {
    
    out <- list(expl = expl)
    
  } else {
    
    delta  <- Fout$delta
    si     <- Fout$si        				  
    s      <- rowMeans(si)  				  

    alpha  <- zeta[2]
    deltar <- delta + alpha*p
    b <- XZPZ %*% deltar                                  
    
    aim         <- -(alpha + zeta[1]*v)
    dsdp        <- -(aim * si) %*% t(si) / NS             
    diag(dsdp)  <- rowMeans(aim * si * (1 - si))          
    
    TH          <- H * dsdp                               
    mkup        <- -as.numeric(solve(TH, tol = 1e-40) %*% s)
    cm          <- p - mkup                               
    
    g <- WZPZ %*% cm         				  
    
    xi <- deltar - XPF %*% b 				  
    wo <- cm - WPF %*% g     				  
    
    out <- list(b = b, g = g, delta = delta, cm = cm, si = si, xi = xi, wo = wo, expl = expl)
  }
  
  return(out)
}

gmm <- function(zeta) {
  
  print(zeta) 			   
  
  strc <- struct(zeta)
  expl <- strc$expl
  
  if(expl == 'try-error') {
    out <- 1e+10 		   
  } else {
    
    xi <- strc$xi
    cm <- strc$cm
    wo <- strc$wo

    wo <- wo + ((cm > p)*(cm - p)^2 + (cm < 0)*(100*cm)^2) 
    
    om  <- c(xi, wo) 		   
    out <- t(om) %*% ZPZ2 %*% om   
    out <- as.numeric(out)/1e+5    
  }
  
  print(out)
  return(out)
}


omega1 <- function(beta) {
  alpha   <- zeta[2]
  delta_r <- delta + alpha*p
  xi <- delta_r - XPF %*% beta
  wo <- wo
  return(c(xi, wo))
}

omega2 <- function(zeta) {
  out  <- struct(zeta)
  xi   <- out$xi
  wo   <- out$wo
  return(c(xi, wo))
}

mktsimul <- function(p, p0, deltar, cm, H, v, zeta, pup) {
 
  alpha <- zeta[2]
  aim   <- -(alpha + zeta[1]*v)  

  mui   <- - zeta[1]*v*p
  delta <- deltar - alpha*p


  s <- share(delta = delta, mui = mui)
  si  <- s$si 
  s   <- s$s

  dsdp        <- -(aim * si) %*% t(si) / NS        
  diag(dsdp)  <- rowMeans(aim * si * (1 - si))      
  TH          <- H * dsdp
  mkup        <- -as.numeric(solve(TH, tol = 1e-40) %*% s)
  eq          <- p - cm - mkup 

  return(eq)
}

p_simul <- function(p, pup){
   p <- pmin(p, (1 + pup) * p0)
   p <- pmax(p, (1 - pup) * p0)
   p <- pmax(p, 10)
   return(p)
}


best_response_2 <- function(p_optim) {
   iter_count <<- iter_count + 1
   p <- numeric(nrow(a))
   p[idx_change] <- p_optim
   p[idx_fixed]  <- p_fixed

   mui    <- - zeta[1]*v*p
   deltar <- a$deltar
   delta  <- deltar - zeta[2]*p
   
   s <- share(delta = delta, mui = mui)
   si  <- s$si 
   s   <- s$s

   s_optim <- s[idx_change]
   M <- a$M[idx_change]
   cm <- a$CM[idx_change]

   pi_optim <- (p_optim - cm) %*% (M*s_optim)
   cat(sprintf("Iter: %d | Beneficio: %.4f | Precio Promedio: %.2f\n", 
               iter_count, pi_optim, mean(p_optim)))
   if(!is.finite(pi_optim)) return(-1e10)
   return(as.numeric(pi_optim)) 
}

welfare <- function(p, deltar, v, zeta, M){
	aim <- zeta[2] + zeta[1]*v
	aim <- aim[1:length(ts), ] 

	mui <- -zeta[1]*v*p
	delta  <- deltar - zeta[2]*p
	
	num <- exp(delta + mui) 
	x   <- t(sapply(unique(t), function(m) colSums(num[which(t == m), ])))
	
	W  <- sum(M*rowMeans(log(x)/aim, na.rm = T))/1000
	return(W)
}

